#include <assert.h>
#include <stdio.h>

#include "../src/mstep.c"

typedef struct _MstepTestDatum MstepTestDatum ;
struct _MstepTestDatum 
{ 
  struct
    {
      float     kappa;
      float     logalpha;
      float     tau[(5-1)*4];
      float     loglambda[4];
      float     pz[5];
      Rating    ratings[8];
      float     logalphastddev;
    }                   inputs;

  struct
    {
      float     newq;
      float     logalpha;
    }                   desired_output;
};

static void
test_mstep (void)
{
  MstepTestDatum data[] = 
    {{{0.405252,  0.9456, {0.73181, 0.750142, 0.745328, 0.567453, 0.574857, 0.502376,  0.755661, 0.533463, 0.0317735, 0.66278, 0.431914, 0.890441,  0.351659, 0.511797, 0.85954, 0.230514}, {0.963732, 0.297182,  0.265408, 0.440542}, {0.160673, 0.486395, 0.860156, 0.494942,  0.428863}, {{0, 0}, {0, 3}, {1, 1}, {1, 3}, {2, 3}, {2, 2}, {2,  4}, {3, 2}},  0.736253}, {-36.341188812110168530, -1.7037836474824718488}},  {{0.114829,  0.927489, {0.854005, 0.233877, 0.359168, 0.394026, 0.822232,  0.571097, 0.927253, 0.503584, 0.470573, 0.0593, 0.0677133,  0.273071, 0.506841, 0.762118, 0.802305, 0.832528}, {0.346169,  0.275723, 0.942149, 0.337586}, {0.917306, 0.53947, 0.82732,  0.410097,  0.0633005}, {{0, 1}, {2, 4}, {0, 3}, {1, 0}, {0, 4}, {2, 2}, {2,  1}, {3, 4}},  0.305593}, {-60.213688728625251074, -1.4850661886454523106}},  {{0.468153,  0.0160714, {0.241069, 0.734496, 0.540899, 0.512487, 0.770495,  0.675196, 0.473186, 0.239416, 0.263654, 0.913078, 0.670881,  0.406888, 0.917486, 0.637354, 0.728732, 0.0693021}, {0.000179541,  0.0978844, 0.901412, 0.659205}, {0.936879, 0.792291, 0.433259,  0.643134,  0.69581}, {{0, 0}, {1, 4}, {0, 4}, {2, 3}, {2, 0}, {1, 1}, {1,  0}, {2, 4}}, 0.0577954}, {-180.36304620840527413,  0.034721291249970015721}}, {{0.89236,  0.130647, {0.925315, 0.3826, 0.419174, 0.89123, 0.661661, 0.469522,  0.748294, 0.484342, 0.744175, 0.832168, 0.0195621, 0.41504,  0.743996, 0.734283, 0.11815, 0.755835}, {0.807117, 0.941992,  0.684891, 0.112701}, {0.111306, 0.884196, 0.792531, 0.982055,  0.185991}, {{3, 2}, {1, 2}, {0, 3}, {0, 4}, {3, 1}, {1, 2}, {2,  4}, {1, 1}},  0.501597}, {-53.363566507665113232, -1.4447334394755153844}},  {{0.373356,  0.0908247, {0.52433, 0.0320748, 0.625062, 0.606483, 0.780155,  0.199907, 0.6055, 0.191443, 0.0361593, 0.465624, 0.48735,  0.435608, 0.229043, 0.523632, 0.802459, 0.322906}, {0.117736,  0.639436, 0.00992804, 0.340851}, {0.931745, 0.137839, 0.636572,  0.250027,  0.407414}, {{1, 1}, {1, 3}, {0, 3}, {1, 0}, {3, 4}, {0, 3}, {3,  1}, {0, 0}},  0.105764}, {-77.564458190383898596, -0.17295969059314759205}},  {{0.0115093,  0.643544, {0.627259, 0.905857, 0.406009, 0.452101, 0.5911,  0.440233, 0.918659, 0.0164935, 0.362057, 0.916601, 0.1162,  0.693587, 0.244321, 0.277165, 0.106272, 0.352736}, {0.312576,  0.139326, 0.469701, 0.102709}, {0.905162, 0.0335611, 0.458191,  0.459165,  0.277903}, {{2, 3}, {3, 3}, {3, 4}, {0, 1}, {2, 2}, {0, 1}, {2,  2}, {2, 1}},  0.127704}, {-45.971365670470138661, -0.39695708516286353640}},  {{0.0521821,  0.00706384, {0.686803, 0.687471, 0.133523, 0.99057, 0.324745,  0.77087, 0.0173226, 0.296983, 0.0804241, 0.493705, 0.91105,  0.944247, 0.767848, 0.35438, 0.44135, 0.841538}, {0.862686,  0.320819, 0.983159, 0.382373}, {0.584783, 0.193115, 0.930976,  0.375309,  0.89798}, {{1, 2}, {2, 3}, {0, 1}, {3, 0}, {3, 3}, {3, 0}, {2,  3}, {2, 0}},  0.505644}, {-48.943220760112974149, -1.8790430532466969382}},  {{0.797453,  0.384739, {0.573235, 0.734774, 0.780131, 0.0877559, 0.492811,  0.241069, 0.869081, 0.143508, 0.724963, 0.886689, 0.427731,  0.30197, 0.862277, 0.56587, 0.444572, 0.919597}, {0.277494,  0.372756, 0.513596, 0.544287}, {0.379514, 0.867112, 0.716142,  0.159548,  0.806279}, {{2, 0}, {2, 4}, {2, 4}, {0, 2}, {0, 3}, {2, 1}, {1,  1}, {3, 4}},  0.132338}, {-118.79300571825089910, -0.29655878979308960848}},  {{0.936012,  0.0717923, {0.313468, 0.89127, 0.066931, 0.928284, 0.588505,  0.00458081, 0.6392, 0.626314, 0.726228, 0.43871, 0.194628,  0.706717, 0.448733, 0.0659547, 0.681032, 0.16243}, {0.0692195,  0.198843, 0.96489, 0.00288165}, {0.262941, 0.0665046, 0.0288779,  0.931089,  0.949473}, {{3, 1}, {1, 1}, {3, 1}, {2, 0}, {2, 4}, {2, 2}, {2,  4}, {3, 4}},  0.175235}, {-89.005668494572267811, -0.38617449859362787234}},  {{0.961947,  0.00280539, {0.360968, 0.170654, 0.322747, 0.376492, 0.63474,  0.731944, 0.128119, 0.669774, 0.186007, 0.665989, 0.447087,  0.507345, 0.116787, 0.467146, 0.482197, 0.504463}, {0.853847,  0.400642, 0.453319, 0.573374}, {0.904374, 0.225407, 0.491372,  0.570568,  0.543406}, {{3, 0}, {0, 4}, {2, 1}, {3, 4}, {2, 3}, {0, 3}, {0,  4}, {3, 1}}, 0.0547526}, {-236.37571494163540496,  0.42877091943914885519}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_logalpha = 
        parameterized_normal_distribution (data[i].inputs.logalphastddev);

      float hq = 0;
      for (unsigned int k = 0; k < 5; ++k)
        {
          hq -= data[i].inputs.pz[k] * log (data[i].inputs.pz[k]);
        }

      float logalpha = data[i].inputs.logalpha;
      float q = mstep (data[i].inputs.kappa,
                       &logalpha,
                       5,
                       data[i].inputs.tau,
                       data[i].inputs.loglambda,
                       data[i].inputs.pz,
                       data[i].inputs.ratings,
                       8,
                       &prior_logalpha.base,
                       hq,
                       1000,
                       1e-7);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-3 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%g ?= %g\n", q, data[i].desired_output.newq),
                0 ));
      
      assert (fabs (logalpha - data[i].desired_output.logalpha) <=
              1e-3 * (1 + fabs (logalpha) + fabs (data[i].desired_output.logalpha)) ||
              ( fprintf (stderr, "(oops ... but this is ok) %g ?= %g\n", logalpha, data[i].desired_output.logalpha),
                1 ));

    }
}

int 
main (void)
{
  srand48 (69);

  test_mstep ();

  return 0;
}
